\documentclass{ctexart} %使用 ctex 字体
\usepackage{xltxtra} % 让 xelatex 支持UTF
\usepackage{graphicx} % 图形库
\usepackage{amsmath} % AMS 的数学符号
\usepackage{listings}


\begin{document}
\title{项目作业}


\author{聂宇潇 \\ 数学与应用数学 3220101556}
\maketitle

\section{有限元方法的基本设置}
这是我们实际使用有限元方法计算某物的第一个示例。我们将解决一个简化版的泊松方程，其中边界值为零，但右手边有非零值：
\begin{align}
-\Delta u &= f\ \ \ \ in\ \Omega, \nonumber \\
u &= 0\ \ \ \ on\ \partial \Omega. \nonumber
\end{align}\notag
我们将在单位正方形$\Omega=[0,1]^2$上解决这个方程.在这个程序中，我们只考虑特殊情况$f(x)=1$。\\
如果你已经学习了有限元方法的基础知识，你会记得我们需要采取哪些步骤来用有限维近似来逼近解u。具体来说，我们首先需要推导上述方程的弱形式，这是通过将方程从左边乘以一个测试函数$\varphi$得到的（我们将在下面解释为什么从左边而不是从右边相乘），然后在域$\Omega$上进行积分：
\[-\int_{\Omega}\varphi\Delta u = \int_{\Omega}\varphi f. \]
这可以通过部分积分来进行积分：
$$\int_{\Omega}\nabla\varphi\cdot\nabla u - \int_{\partial\Omega}\varphi n\cdot\nabla u = \int_{\Omega}\varphi f.$$
测试函数$\varphi$必须满足相同类型的边界条件（数学上的术语：它必须来自于我们寻求解的集合的切空间），因此在边界上$ \varphi=0$，因此我们要找的弱形式可以表示为：
$$(\nabla\varphi,\nabla u) = (\varphi,f),$$
在这里，我们使用了常见的记号 $(a,b)=\int_{\Omega}ab$。然后，问题要求找到一个函数 u，对于所有来自适当空间（这里是 $H^1$空间）的测试函数 $\varphi$，该表达式都成立。\\
当然，在一般情况下我们无法在计算机上找到这样的函数，相反我们寻求一个近似解 $u_h(x)=\sum_jU_j\varphi_j(X)$，其中$ U_j $是需要确定的未知展开系数（该问题的“自由度”），而$ \varphi_i(x)$ 是我们将要使用的有限元形状函数。为了定义这些形状函数，我们需要以下信息：
\begin{itemize}
  \item 一个用于定义形状函数的网格。你已经在步骤1和步骤2中看到了如何生成和操作描述网格的对象。
  \item 一个有限元，描述了我们要在参考单元上使用的形状函数（在deal.II中，参考单元总是单位区间 $[0,1]$、单位正方形 $[0,1]^2$ 或单位立方体$ [0,1]^3$，具体取决于你所工作的空间维度）。在步骤2中，我们已经使用过类型为 $FE\_Q<2>$ 的对象，它表示通常的拉格朗日元素，通过在支持点上进行插值来定义形状函数。最简单的是 $FE\_Q<2>(1)$，它使用一次多项式。在二维中，它们通常被称为双线性，因为它们在参考单元的每个坐标上都是线性的。（在一维中，它们将是线性的，在三维中将是三线性的；然而，在 deal.II 的文档中，我们通常不做这个区分，简单地将这些函数都称为“线性”。)
  \item 一个 DoFHandler 对象，枚举网格上的所有自由度，以有限元对象提供的参考单元描述为基础。你在步骤2中已经看到了如何实现这一点。
  \item 一个映射，用于说明如何从有限元类在参考单元上定义的形状函数获得实际单元上的形状函数。默认情况下，除非你明确指定，deal.II将使用一个（双线性、三线性）映射进行此操作，所以在大多数情况下，你不必担心这一步骤。
\end{itemize}
通过这些步骤，我们现在有一组函数 $\varphi_i$，并且我们可以定义离散问题的弱形式：找到一个函数 $u_h$，即找到上面提到的展开系数 $U_j$，使得
\[(\nabla \varphi_i,\nabla u_h) = (\varphi_i,f),i=0...     N-1.\]
请注意，我们这里遵循从零开始计数的惯例，这在 C 和 C++ 中很常见。如果将表示 $u_h(x) = \sum_jU_j\varphi_j(x)$ 插入到方程中，那么这个方程可以重写为一个线性系统，并且观察到
\begin{align}
(\nabla\varphi_i, \nabla u_h) &= (\nabla\varphi_i, \nabla[\sum_j U_j\phi_j]) \\
&= \sum_j(\nabla\varphi_i, \nabla[U_j\varphi_j]) \\
&= \sum_j(\nabla\varphi_i, \nabla\varphi_j)U_j
\end{align}\notag
有了这个，问题可以表述为：找到一个向量 U，使得\[AU = F,\]
其中矩阵 A 和右边界 F 的定义如下:
\begin{align}
  A_{ij} &= (\nabla\varphi_i,\nabla\varphi_j)\\
  F_i &= (\varphi_i,f)
\end{align}\notag

\section{我们应该左乘还是右乘测试函数？}
在我们继续描述如何计算这些量之前，请注意，如果我们将原始方程从右边而不是左边乘以一个测试函数，那么我们将得到以下形式的线性系统：
$$U^TA = F^T$$
通过转置$F^T$，这当然等价于解决以下形式的问题：
$$A^TU = F$$
这在这里与上面是相同的，因为 $A = A^T$。但一般情况下并非如此，为了避免任何混淆，经验表明，习惯上从左边而不是从右边乘以方程（通常在数学文献中这样做）可以避免一类常见的错误，因为矩阵在比较理论和实现时会自动正确，无需进行转置。请参考本教程中的第一个示例步骤-9，我们在其中有一个非对称的双线性形式，从右边或从左边乘以方程会产生不同的结果。

\section{计算矩阵和右手边向量的方法}
现在我们知道我们需要什么（即：保存矩阵和向量的对象，以及计算$ A_{ij}$、$F_i$ 的方法），我们可以看看实现这一点需要做些什么：
\begin{itemize}
  \item A、U和F的对象类型为SparseMatrix和Vector，我们将在下面的程序中看到用于解决线性系统的类。
  \item 我们需要一种形成积分的方法。在有限元方法中，最常用的方法是使用数值积分，即将积分替换为对每个单元上一组点的加权求和。也就是说，我们首先将关于$\Omega$的积分分解为所有单元上的积分：
    \begin{align}
      A_{ij} &= (\nabla\varphi_i,\nabla\varphi_j) = \sum_{K\in T}\int_K\nabla\varphi_i\cdot\nabla\varphi_j,\\
      F_i &= (\varphi_i,f) = \int_{K\in T}\int_K\varphi_if,
    \end{align}\notag
    然后通过数值积分来近似每个单元的贡献：
    \begin{align}
      A_{ij}^K &= \int_K\nabla\varphi_i\cdot\nabla\varphi_j \approx \sum_q\nabla\varphi_i(x_q^K)\cdot\nabla\varphi_j(x_q^K)\omega_q^K,\\
      F_i^K &= \int_K\varphi_if \approx \sum_q\varphi_i(x_q^K)f(x_q^K)\omega_q^K,
    \end{align}\notag
    其中，$x_q^K$表示单元K上的第 q 个数值积分点，$\omega_q^K$表示第 q 个数值积分权重。这个过程需要涉及到多个方面，下面我们将逐一讨论它们。
  \item 首先，我们需要一种方法来描述数值积分点$x_q^K$的位置和对应的权重$\omega_q^K$。通常，它们会通过与形函数相同的方式从参考单元映射得到，也就是隐式地使用MappingQ1类，或者如果你明确指定，通过其他派生自Mapping的类进行映射。参考单元上的积分点和权重由继承自Quadrature基类的对象来描述。通常，我们可以选择一个数值积分公式（即一组点和权重），使得数值积分能够精确等于矩阵的积分；这可以通过高斯数值积分公式来实现，它在QGauss类中实现了。
  \item 然后，我们需要一些工具来帮助我们在单元K上评估$\varphi_i(x_q^K)$。这就是FEValues类的作用：它接受一个有限元对象来描述参考单元上的$\varphi$，一个数值积分对象来描述数值积分点和权重，以及一个映射对象（或者隐式使用MappingQ1类），并在实际单元K上提供形函数在数值积分点上的值和导数，以及其他用于积分的各种信息。
\end{itemize}
FEValues在组装过程中真正起到了核心作用。一个可以理解FEValues的方式是：FiniteElement及其派生类描述了形函数，也就是无限维对象：函数在每个点都有值。我们出于理论上的需要，希望使用函数的积分来进行分析。然而，对于计算机来说，这是一个非常困难的概念，因为它们通常只能处理有限数量的信息，所以我们通过映射（映射对象）和参考单元上定义的点（数值积分对象），将积分替换为数值积分点的求和，将参考单元上的点映射到实际单元上的点。本质上，我们将问题简化为只需要有限数量的信息，即形函数值和导数、数值积分权重、法向量等，仅在有限的一组点上。FEValues类将这三个组件结合在一起，并在特定单元K上提供了这个有限集合的信息。当我们在下面组装线性系统时，你将看到它发挥的作用。\\
值得注意的是，如果您在应用程序中自己创建这三个对象，并自己处理信息，也可以实现所有这些。然而，这既不会更简单（FEValues类提供了您实际需要的精确信息），也不会更快：FEValues类经过高度优化，只计算您需要的特定单元上的信息；如果任何东西可以从先前的单元中复用，它将这样做，并且该类中有很多代码来确保在有利的情况下进行缓存。\\
介绍的最后一部分是提到，在得到线性系统之后，我们使用迭代求解器来解决它，然后进行后处理：我们使用DataOut类创建一个输出文件，然后可以使用常见的可视化程序对其进行可视化。
\section{关于实现的部分}
虽然这是使用有限元方法求解的最简单的方程，但该程序展示了大多数有限元程序的基本结构，并且也作为接下来几个程序的基本模板。具体而言，该程序的主要类如下所示：
\begin{lstlisting}[language=C++, caption=class Step3]
  class Step3
{
  public:
    Step3 ();
    void run ();
  private:
    void make_grid ();
    void setup_system ();
    void assemble_system ();
    void solve ();
    void output_results () const;
    Triangulation<2>     triangulation;
    FE_Q<2>              fe;
    DoFHandler<2>        dof_handler;
    SparsityPattern      sparsity_pattern;
    SparseMatrix<double> system_matrix;
    Vector<double>       solution;
    Vector<double>       system_rhs;
};
\end{lstlisting}
这遵循了面向对象编程的数据封装原则，即我们尽力将该类的大部分内部细节隐藏在对外不可访问的私有成员中。\\
让我们从成员变量开始：这些成员变量遵循我们在上面的要点中概述的构建模块，即我们需要一个Triangulation和一个DoFHandler对象，以及一个描述我们想要使用的形状函数类型的有限元素对象。第二组对象与线性代数相关：系统矩阵、右手边以及解向量，以及描述矩阵的稀疏模式的对象。这是该类所需的全部内容（也是任何用于静态PDE求解器所需的基本要素），并且需要在整个程序中保留下来。与此相反，我们在装配过程中仅需要FEValues对象，并且将其作为局部对象在执行装配的函数中创建，然后在函数结束时销毁。\\
其次，让我们来看一下成员函数。同样，这些函数已经形成了几乎所有后续教程程序将使用的常见结构：
\begin{itemize}
  \item $make\_grid()$:这可以称为预处理函数。顾名思义，它设置了存储三角剖分的对象。在后续示例中，它还可以处理边界条件、几何形状等内容。\\
  \item $setup\_system()$:然后，在这个函数中设置了解决问题所需的所有其他数据结构。特别地，它将初始化DoFHandler对象并正确调整与线性代数相关的各种对象的大小。这个函数通常与上面的预处理函数分开，因为在时间相关的程序中，每当网格被自适应细化(我们将在step-6中看到)时，它可能至少在每几个时间步骤中被调用一次。另一方面，在上面的预处理函数中设置网格本身仅在程序开始时执行一次，因此被分离成独立的函数。
  \item $assemble\_system()$:因此，这是在上面的介绍中详细讨论过的计算矩阵和右手边的内容的地方。由于对这个线性系统进行操作在概念上与计算其条目非常不同，我们将其与后续的函数分开。
  \item $solve()$:因此，这个函数是用来计算线性系统$AU=F$的解U的函数。在当前程序中，由于矩阵非常简单，这是一个简单的任务。但是，当问题变得不再那么简单时（例如，在了解更多关于该库的内容后，可以参考step-20、step-22或step-31），计算解U将成为程序规模的重要组成部分。
  \item $output\_results()$:最后，当你计算得到一个解时，你可能希望对它进行一些处理。例如，你可能想以可视化的格式输出它，或者你可能想计算你感兴趣的量：比如热交换器中的热流量、翼的空气摩擦系数、桥梁的最大荷载，或者仅仅是数值解在某一点上的值。因此，这个函数是用来进行后处理的地方。
\end{itemize}
所有这些都由唯一的公共函数（除了构造函数）连接在一起，即run()函数。这个函数是从创建此类型对象的地方调用的，它按照正确的顺序调用所有其他函数。将这个操作封装到run()函数中，而不是从main()函数调用所有其他函数，可以确保你可以改变这个类内部的关注点分离是如何实现的。例如，如果其中一个函数变得太大，你可以将它分割成两个函数，你需要关注更改的只有这个类内部的地方，而不是任何其他地方。\\
正如上面提到的，你将在许多接下来的教程程序中再次看到这个通用的结构，有时函数的名称拼写可能会有所变化，但基本上是按照这个功能分离的顺序。

\section{关于类型的注释}
在命名空间types中，deal.II通过别名定义了一些整数类型。特别是，在这个程序中你会在几个地方看到$types::global\_dof\_index$这是一个整数类型，用于表示自由度的全局索引，即在一个基于三角剖分的DoFHandler对象中表示特定自由度的索引（与在特定单元内表示特定自由度的索引相对）。对于当前程序（以及几乎所有的教程程序），你将有几千到几百万个全局未知数（对于$Q_1$元素，在2D中每个单元上有4个本地未知数，在3D中有8个）。因此，一个足够大的数据类型来存储全局自由度索引的数字是unsigned int，因为它可以存储0到略大于40亿的数字（在大多数系统上，整数是32位）。事实上，这就是$types::global\_dof\_index$的定义。\\
那么，为什么不直接使用unsigned int呢？在deal.II的7.3版本之前，它确实是这样做的。然而，deal.II支持非常大规模的计算（通过第40步中讨论的框架），当这些计算分布在几千个处理器上时，可能会涉及到超过40亿个未知数。因此，在某些情况下，unsigned int并不足够大，我们需要一个64位的无符号整数类型。为了实现这一点，我们引入了$types::global\_dof\_index$，默认情况下它被定义为unsigned int，但如果有必要，可以通过在配置过程中传递特定的标志来将其定义为unsigned long long int。\\
这涵盖了技术方面的内容。但还有一个文档目的：在库和基于它构建的代码中，如果你看到一个使用数据类型$types::global\_dof\_index$的地方，你就立即知道被引用的是一个全局自由度索引。如果我们只是使用unsigned int（它也可能是一个局部索引、边界指示器或材料标识等），就不会有这样的含义明显。立即知道变量所指的是什么也有助于避免错误：如果你看到一个类型为$types::global\_dof\_index$的对象被赋值给类型为$types::subdomain\_id$的变量，尽管它们都由无符号整数表示，编译器不会抱怨，但很明显肯定存在一个错误。\\
更具实际意义的是，在装配过程中，我们会创建一个4×4的矩阵（在2D中使用$Q_1$单元），用来表示当前单元格的贡献，然后我们需要将这个矩阵的元素添加到全局（系统）矩阵的相应元素中。为此，我们需要获取与当前单元格相关的局部自由度的全局索引，对于这一点，我们将始终使用以下代码片段：
\begin{lstlisting}[language=C++, caption=]
  cell->get\_dof\_indices (local\_dof\_indices);
\end{lstlisting}
其中，$local\_dof\_indices$被声明为：
\begin{lstlisting}[language=C++, caption=]
 std::vector<types::global\_dof\_index> local\_dof\_in
 dices (fe.dofs\_per\_cell);
\end{lstlisting}

这个变量的名称可能有点不准确，它代表着“在当前单元格上局部定义的自由度的全局索引”，但在整个库中，保存这些信息的变量都普遍以这种方式命名。

\section{step-3所计算的问题及所用参数}
Step-3介绍了一个简单的Poisson方程求解器。该教程的目标是求解具有Dirichlet边界条件的二维Poisson方程。以下是问题和具体参数的描述：\\
问题描述：\\
在二维域上求解poisson方程
\[-\Delta u = f \]
其中u是未知函数，f是已知函数。\\
具体参数：
\begin{itemize}
\item 域的大小和几何形状：通过设置域的大小和形状来定义求解区域的几何结构。\\
\item 网格剖分：用于离散化求解区域的有限元网格，可以设置不同的细粒度、网格大小等参数来控制网格的精细程度。\\
\item 有限元空间：选择合适的有限元空间来近似解，例如选择连续Galerkin有限元空间。\\
\item 弱形式：构建求解离散化问题的弱形式，例如使用变分原理得到离散化问题的等效形式。\\
\item 边界条件：规定边界上变量的值或导数的值。在Step-3中，以线性插值的方式在边界上规定 u 的值。\\
\item 数值求解器：选择合适的求解方法和参数，例如使用共轭梯度法或GMRES方法等\\
\end{itemize}

\section{运行step-3产生数据的图像}
运行step-3，将所得数据利用gnulpot进行绘图，可得到下图：
\begin{figure}
  \centering
  \includegraphics[width=0.5\textwidth]{tu.pdf}
  \caption{step-3}
  \label{fig:1}
\end{figure}






\end{document}
